program newton_polynomial
  implicit none
  real,allocatable,dimension(:)::x,f,a
  real,allocatable,dimension(:,:)::D
  integer::i,j,n
  real::xp,fp,su,pi

  print*,'enter number of data points'
  read*,n
  allocate(x(n))
  allocate(f(n))
  allocate(a(n))
  allocate(D(n,n))
  print*, 'enter co-ordinates and function values(x(i) y(i))'
  do i = 1,n
    read*,x(i), f(i)
  end do
  !constructing a difference table D

  do i = 1,n
    D(i,1) = F(i)
  end do

  do j = 2,n
    do i =1,n-j+1
      d(i,j) = (D(i+1,j-1) - D(i,j-1))/(x(i+j-1)-x(i))
    end do
  end do

  !set the coefficients of interpolating polynomial
  do j = 1,n
    a(j) = D(i,j)
  end do

  !compute interpolation values
  print*,'enter co-ordinate'
  read*,xp
  su = a(1)
  do i = 2,n
    pi = 1.0
    do j = 1, i-1
      pi = pi*(xp - x(j))
    end do
    su = su + a(i)*pi
  end do
  fp = su

  print*,'The interpolated value is:',fp
end program




